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This article describes an error covariance analysis methodology used to investi- 
gate different weighting schemes for two-way ( coherent ) Doppler data in the presence 
of transmission-media and observing-platform calibration errors. The analysis fo- 
cuses on orbit-determination performance in the interplanetary cruise phase of deep- 
space missions. Analytical models for the Doppler observable and for transmission- 
media and observing-platform calibration errors are presented, drawn primarily 
from previous work. Previously published analytical models have been improved 
upon by (1) considering the effects of errors in the calibration of radio signal prop- 
agation through the troposphere and ionosphere as well as station-location errors; 

(2) modelling the spacecraft state transition matrix using a more accurate piecewise- 
linear approximation to represent the evolution of the spacecraft trajectory; and 

(3) incorporating Doppler data weighting functions that are functions of elevation 
angle, which reduce the sensitivity of the estimated spacecraft trajectory to tropo- 
sphere and ionosphere calibration errors. The analysis is motivated by the need to 
develop suitable weighting functions for two-way Doppler data acquired at 8.4 GHz 
(X-band) and 32 GHz (Ka-band). This weighting is likely to be different from that 
in the weighting functions currently in use; the current functions were constructed 
originally for use with 2.3-GHz (S-band) Doppler data, which are affected much 
more strongly by the ionosphere than are the higher frequency data. 


I. Introduction 

Interplanetary spacecraft navigation is accomplished 
using a number of different techniques. The classical 
means of navigation have included those based on measure- 
ments of Doppler shift and on ranging measurements. In 


recent years, differential very long baseline interferometry 
(VLBI) has been added to the arsenal of tools used for nav- 
igation [1]. For planetary approach and encounter, optical 
navigation using images of solar-system bodies against the 
background stars also is a very useful means of naviga- 
tion [2]. Differenced (two-way minus three-way) Doppler 
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and range currently are being investigated for spacecraft 
tracking [3,4], as are other interferometric data types such 
as connected-element interferometry [5] and same-beam 
VLBI [6]. 

Despite the plethora of data types being developed, the 
workhorse navigation data type remains two-way (coher- 
ent) Doppler. These data are easy to acquire and can 
be obtained without interference with telemetry transmis- 
sion from a spacecraft. However, fairly long tracking ses- 
sions (“passes”) are necessary on each day of tracking. 
Doppler data are directly sensitive to motion along the 
line of sight; accumulation of Doppler data over a num- 
ber of days provides increasing sensitivity to position and 
velocity in the plane of the sky (the plane perpendicular 
to the Earth-spacecraft line-of-sight). VLBI data yield- 
ing similar or improved angular accuracies often can be 
acquired in far less time at a greater operational cost and 
at the expense of interfering with telemetry transmission 
from the spacecraft. As instrumentation noise has been 
reduced and the calibration of various effects that cause 
systematic measurement errors has advanced, the orbit- 
determination performance obtained from Doppler data 
has improved substantially. 

Obtaining the best accuracy from Doppler data often 
has required that an empirical trade-off be made between 
acquiring the most data possible in a given pass and down- 
weighting or deleting the data which can contain large er- 
rors due to imperfect calibration. In particular, it has been 
a common practice to deweight or discard data acquired 
at low elevation angles because of errors in the calibration 
of propagation effects in the Earth’s ionosphere and tro- 
posphere. For Doppler data taken near 2.3 GHz (S-band), 
ionosphere calibration errors typically dominate the low- 
elevation data. Current missions often use the higher radio 
frequency of 8.4 GHz (X-band), and future missions may 
use transmission frequencies near 32 GHz (Ka-band) for 
improved telemetry performance. At these higher frequen- 
cies, the ionosphere calibration errors become less impor- 
tant (ionospheric propagation errors scale with the inverse 
square of the radio frequency), and troposphere calibra- 
tion uncertainties dominate the errors in the low-elevation 
data. For Doppler data used in orbit determination, the 
weighting function currently implemented at JPL was es- 
tablished when radio metric data were obtained at 2.3 GHz 
and emphasized the minimization of errors caused by im- 
perfect ionosphere calibration. Therefore this function is 
inappropriate for 8.4-GHz and 32-GHz data. 

This article describes the development of a numerical 
error covariance analysis used to investigate the benefits 
of new elevation-dependent weighting schemes for Doppler 


data acquired at high radio frequencies. The analysis de- 
scribed below is applicable to interplanetary cruise trajec- 
tories, but does not include features to enable computation 
of gravitational effects encountered during planetary ap- 
proach or encounter phases. A piecewise-linear approxima- 
tion to the evolution of the spacecraft state vector enables 
the analysis to be applied to data arcs that are months in 
length; keeping quadratic terms could enable the study of 
even longer data arcs at the expense of much additional 
computation. The virtue of the piecewise-linear analysis 
is that it is relatively simple, enabling consideration of a 
large number of different cases in a short period of time. 

Previous analytical and numerical studies of Doppler- 
based orbit-determination performance can be found in 
a number of different references. For example, Hamilton 
and Melbourne [7] and Thurman [8] have considered the 
information contained in a single pass of Doppler mea- 
surements. Analyses have been extended to several days 
of data in a number of reports (e.g., [9-11]). However, 
most of these analyses have considered only data noise 
and, in some cases, station-location errors. The analysis 
and computational procedures described here are based 
on that previous work, but incorporate a number of im- 
provements. A better approximation hats been used for 
the propagation of observables from observation times to 
the reference epoch. Some terms that have been neglected 
in the past have been included. Effects of the errors in 
calibration of the static troposphere, the ionosphere, and 
the station locations all have been considered simultane- 
ously, with special attention paid to the accuracy of the 
troposphere and ionosphere models. Finally, the method- 
ology includes the ability to use different Doppler data- 
weighting functions in order to enable investigation of the 
improvement in orbit-determination accuracy that might 
be achieved by using different weighting schemes. 

II. Summary of Least-Squares Error 
Covariance Analysis 

The description of least-squares covariance analysis is 
kept to a minimum here; the reader should consult other 
references (e.g., [7,8,12]) for more explanation. Consider 
the range-rate pi associated with the tth Doppler mea- 
surement. This measurement is modeled as a function of 
the epoch trajectory ro, the measurement time and the 
random measurement error n,-, as follows: 

Pi = f(ro,U ) + (1) 

In a linearized model, measurement residuals A p, can be 
related to perturbations Afo in the trajectory by 
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A pi 


dpi _ 

—r A r 0 + rii 

dr 0 


( 2 ) 


The accumulated information matrix, J, for a set of N 
independent measurements of pi with individual weights 
Wi is 


N 

J = ^2 Wi ( dpi/drof {dpi/dro) (3) 

»=i 

If each data point were weighted equally according to an 
estimate of the variance ao 2 of the Doppler measurement, 
Wi would be replaced by cp -2 , which could be taken out- 
side the summation. However, since a primary goal of this 
investigation is to study the effect of a variety of data- 
weighting schemes, the more general formulation of Eq. (3) 
will be maintained. 

Equation (3) holds for points within a continuous track- 
ing pass as well as for points acquired over several track- 
ing passes. If it is assumed that no a priori information is 
available, the error covariance matrix A 2 is the inverse of 
the information matrix: 


A 2 = J~ l 


(4) 


III. Partial Derivatives and the Computed 
Doppler Error Covariance 

Consider the case of a spacecraft whose right ascension 
and declination are denoted by (cc, 6), viewed from a track- 
ing station whose cylindrical coordinates with respect to 
the Earth’s center are ( r,,\,,z , ). The station-spacecraft 
tracking geometry is shown in Fig. 1. The range-rate for 
the spacecraft, p, which can be related directly to the mea- 
sured Doppler shift, is given approximately by [10] 


r — z,8 cos 6 + r a (ip — d) cos 6 sin ( ip — a) 


-I- r,6 sin 8 cos (<p — a) 


(5) 


In Eq. (5), r is the spacecraft distance from the center 
of the Earth, ip (Local Sidereal Time) is the angle be- 
tween the station meridian and the vernal equinox, and 
time derivatives are indicated by dots above the variables. 
The quantity (<p — a) is the hour angle of the spacecraft, 
measured west from the local meridian. The partial deriva- 
tives of the range-rate with respect to the six-parameter 
state vector r = ( r,6,a,r,6,d ) are given by 


/ dp/ dr \ 


/ . 0 
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dp/dS 


z s (5sin 8 — r,(ip — a) sin 8 sin (ip — a) + r t 8 cos 6 cos (<p — 
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( 6 ) 


To relate changes in the measured Doppler frequency 
to changes in the spacecraft state vector ro at some ref- 
erence epoch 1 0 , one must find the spacecraft state tran- 
sition matrix, consisting of the partial derivatives of the 
state vector r with respect to the epoch state vector ro = 
(r o,8o,ao,ro,8o,cto). These partial derivatives constitute a 
6x6 matrix A 2 that is a function of the geometry at the 
reference time and the interval that has passed since that 
reference time. Formally, the matrix A 2 is expressed as 
df/dr 0 . The most common method of deriving this ma- 
trix analytically (e.g., [9]) is to linearize all six components 
of the state vector in time t a s follows: 
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The acceleration terms f'o, ^o, and do in Eq. ( 7 ) can be 
expressed as functions of r 0 , thereby enabling a determi- 
nation of all the relevant partial derivatives. The effects 
of gravitational acceleration due to the Sun are included 
in the computation of the partials. The gravitational ef- 
fects of other bodies have been neglected, so the analysis 
is useful only for the interplanetary cruise portion of a tra- 
jectory, not for the planetary approach/encounter phase. 


Formulas for the entries in the matrix A 2 have been 
derived in [8] and [ 10 ]. Occasional sign errors occur in 
the literature; correct expressions for these partials are 
given in the Appendix. The formulation here remains cast 
in terms of the state vector of the spacecraft instead of 
the coefficients (a,b,c t d,e,f) defined in [6] and [10] to be 
functions of the six parameters in the state vector. 


Most previous covariance analyses (e.g., [ 9 ]) employed 
an A 2 matrix whose terms were functions only of the state 
vector ro at the epoch time and the time interval since that 
epoch time. Such a formulation is a good approximation 
for multiple passes to the extent that the spacecraft motion 
fits the linear model over the entire data arc. 


A 2 - B 2 (t) = M j (t)A i ~ l A i j :\---AlA\ 


= m) fl A itl 


— m 
m 


(8) 


m = 2 


The matrix Aj -1 represents the mapping matrix from the 
midpoint of pass j to the midpoint of pass j — 1 . Compu- 
tation of A} -1 is similar to the computation of the matrix 
A 2 . The straightforward method is to use a linearization 
process similar to that of Eqs. ( 7 ): 
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( 9 ) 


A more accurate representation of the partial deriva- 
tives indicated in Eq. ( 3 ) can be obtained by using a 
piecewise-linear approximation. For each pass, an epoch 
time is chosen to be the time of spacecraft transit during 
that pass. A linear formulation similar to Eqs. ( 7 ) then is 
used to relate the state vector at any time during a pass to 
an epoch state vector for the same pass, using the param- 
eter values for the epoch time of the given pass, not those 
for the epoch time of the first pass. Multiplication by the 
mapping matrices relating each day to the previous day, 
with the appropriate epoch state vector used for each one- 
day step in the propagation, provides the mapping matrix 
to the reference epoch. Hence, the continuously curving 
spacecraft trajectory is approximated by a series of (six- 
dimensional) straight lines of different slope rather than 
by a single straight line. 

Mathematically, let the matrix of partials of the state 
vector for pass j at time t with respect to the epoch state 
vector at the midpoint of pass j be denoted by Mj(t). This 
matrix has the same entries as the matrix A 2 (results given 
in the Appendix), except that r*o is replaced by r} o, the 
spacecraft state vector at the transit time for pass j, and 
the time t is referred to the reference time for pass j, not 
to the first pass. To relate measurements at time t during 
pass j to the epoch state vector on day 1, the matrix A 2 
should be replaced as follows: 


In Eqs. ( 9 ), fj,o is the state vector of the spacecraft at the 
transit time for pass j, is the similar state vector for 

the previous pass, and r is the time between the two tran- 
sits. In the simple piecewise-linear approximation, Aj _1 
would have entries similar to those for A 2 (see Appendix), 
except that f’o would be replaced by fj-1,0 and t would be 
replaced by r. A slightly better approximation has been 
used in the analysis described here: The entries of Aj~ 
are evaluated at fj,j-i = 0 . 5 (f)_i,o + >7, 0) instead of at 
fj_i i0 . Thus, the propagation from the transit time for 
pass j to the transit time for pass j — 1 uses the aver- 
age slope of the six-dimensional state vector over the time 
period between these two times rather than the slope at 
either end of the period. 

Figure 2 is a schematic of several methods of propa- 
gating data to the reference time; the curvature of a tra- 
jectory over 10 days (i.e., the acceleration) has been ex- 
aggerated to highlight the differences. If the actual evo- 
lution of a spacecraft state vector is represented by the 
lower quadratic curve D, the propagation of a state vector 
found for day 10 to the reference time on day 1 can be 
represented by three different approximations to the slope 
of that curve. The horizontal line (labelled A) is the re- 
sult of using the slope on day 1 to estimate the state vector 
that would be found for day 1 by propagating back in time 
from the state vector actually observed on day 10; this is 
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roughly equivalent to propagating observations from day 
10 to the reference time on day 1 , using only the epoch 
state vector of day 1 to evaluate the matrix A 2 . The up- 
permost piecewise-linear “trajectory” (labelled B) uses the 
slope on day j—1 to propagate the trajectory from day j to 
day j — 1 ; it is the same (in principle) as evaluating each 
matrix Aj 1 at the point r)_ 10 (namely, the spacecraft 
state vector at the midpoint of pass j — 1 ) and performing 
the calculation indicated by Eq. ( 8 ). Another piecewise- 
linear approximation (labelled C), virtually indistinguish- 
able from the actual quadratic trajectory (labelled D) in 
Fig. 2, is equivalent to evaluating each Aj _1 at fjy_i (de- 
fined above) and performing the multiplications in Eq. ( 8 ). 
Clearly, the piecewise-linear approximation described in 
connection with Eqs. ( 8 ) and (9) is the best of the three 
methods shown, although the degree to which it surpasses 
the other methods depends on the actual curvature of the 
spacecraft trajectory over a particular data arc. 

The partial derivatives of each Doppler measurement 
with respect to the epoch state vector now can be deter- 
mined readily from 


dpi 

dr 0 


% B ’<‘> 


( 10 ) 


where Z ? 2 is evaluated as shown in Eq. ( 8 ) and described in 
the accompanying text. Using Eqs. (3), (4), and (10), the 
error covariance matrix associated with a weighted least- 
squares estimate of r 0 can be calculated. 

It often is convenient to express the epoch state error 
covariance in a Cartesian coordinate frame rather than 
with spherical coordinates. A particularly useful Cartesian 
frame is the plane-of-sky coordinate system defined by the 
spacecraft position at epoch, which is depicted in Fig. 3. 
The position components r c , ra, and r a represent small 
position perturbations in the direction of increasing r, 6, 
and a coordinate directions relative to the nominal epoch 
spacecraft position. To make the conversion, it is necessary 
to compute the matrix A 3 containing the partial deriva- 
tives of the Cartesian coordinates at the reference time, 
( r c,o, ra,o, '’a.o, fio.o), with respect to the spheri- 

cal coordinates (i- 0 ,6o,a 0 ,r 0 ,S 0 ,a 0 ). Reference [11] gives 
those derivatives, but contains several errors; correct val- 
ues for the partials matrix A 3 are given in the Appendix. 


In Eq. (11), A sp h is simply the error covariance matrix for 
the spacecraft spherical coordinates at epoch. 


IV. Description of Consider-State Error 
Covariance Analysis 

A. Brief Review of Theory 

The error covariance calculation outlined in Sections II 
and III could be used to predict the orbit-determination 
error statistics obtained from a series of Doppler measure- 
ments if each measurement were affected only by random, 
uncorrelated errors. However, systematic errors also can 
be present in the Doppler data due to imperfect calibra- 
tion of a number of quantities that often are not solved 
for in the orbit-determination process. The sensitivity of 
the computed trajectory to these error sources can be esti- 
mated by allowing their influence to be “considered” in de- 
termining the error statistics associated with the computed 
trajectory. Obviously, the use of such a procedure results 
in a suboptimal estimate of the spacecraft trajectory, an 
estimate that does not minimize the expected value of the 
mean-squared estimation errors. This subsection reviews 
briefly the procedures for calculating a “consider” error co- 
variance, while the following subsection describes the par- 
ticular set of considered parameters included in the current 
analysis. 

A consider-state covariance analysis primarily involves 
the computation of the partial derivatives of the estimated 
epoch state vector, ro ie , with respect to the consider pa- 
rameters y c . The consider parameters are those quantities 
that are not estimated explicitly in the orbit determina- 
tion. To perform the computation, the partial derivatives 
of the Doppler observable p with respect to the consider 
parameters y c must be found and then used in combina- 
tion with the partials of p with respect to fo to calculate 
the desired quantities. (Note that ro ie is the estimated 
state vector, while r"o is the actual state vector.) Formally, 
define the matrix dro ie /dy c as the sensitivity matrix, de- 
noted as S. A detailed derivation of the sensitivity matrix 
is given by McReynolds [13]; only the results will be sum- 
marized here. The sensitivity matrix can be expressed as 
a function of the computed covariance, A 2 , representing 
the uncertainty in the estimated parameters due to the 
random measurement errors (the “data noise”), and an 
additional matrix as follows: 


It easily can be shown that the plane-of-sky Cartesian 
covariance matrix, A p0 si is given by 

Apos = A3A S ph^3 (11) 


S = dr 0ie /dy c 

= A 2 £ m ( dpi/dr 0 ) T (dpi/dyc) 

*-l 


} (12) 
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Recall that the computed covariance, A 2 , given in 
Eq. (4) is just the inverse of the information matrix. Now 
let the matrix A y represent the error covariance matrix 
for the consider parameters. The matrix A y is assumed 
to be diagonal in this analysis, with entries equal to the 
variances in the knowledge of the individual consider pa- 
rameters. The statistical uncertainty introduced into the and 
estimated parameters (the spacecraft epoch state vector 
elements) by the consider parameters, denoted by the co- 
variance matrix Ai, is given by 


/Wi 

0 • • 

0 > 


0 

U> 2 • • 

0 




0 

(15c) 
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k 0 

0 • • 

w N ) 


Aq = 

I ■ (To 2 


(15d) 


Ai = SA y S T (13) 

Finally, the total, or consider, error covariance for the esti- 
mated parameters, A3, due to both random measurement 
errors and the consider parameters, is 

A 3 = Ai + A 2 (14) 


For W = I Eq. (15a) reduces to A 2 = (A T WA)~ 1 , 

as it should. To form the sensitivity matrix in the general 
case, Eq. (12) still is appropriate, but with the expres- 
sion in Eq. (15a) used for the computed covariance A 2 . 
The formula for the consider error covariance A 3 , given 
in Eq. (14), then is unchanged after making the necessary 
modification to Eq. (12). This general formulation was 
used in the analysis described herein. 


If desired, the error covariance matrix for the spacecraft 
epoch state vector given in Eq. (14) can be expressed read- 
ily in a plane-of-sky Cartesian frame using the mapping 
formula in Eq. (11). 

The formulation given thus far for the consider error 
covariance is valid only if the weights id,- assigned to the 
Doppler data are equal to 1 /<To 2 - If the weights used are 
not the inverse of the data noise variance, then the formu- 
las given for Ai and A 2 in Eqs. (13) and (4) must be mod- 
ified. This scenario is of particular interest here, since the 
following section describes a Doppler weighting technique 
in which the data weights are not equal to the inverse of 
the data noise variance. This weighting scheme is designed 
to reduce the sensitivity of Doppler-based trajectory esti- 
mates to errors in the calibration of the transmission media 
and the station locations. For this more general problem, 
it can be shown 1 that the computed covariance, A 2 , is 

A 2 = ( A T WA)~ 1 A T WA r ,WA(A T WA ) _1 (15a) 

where 


/ dpi/dfo \ 

I dp 2 /df 0 I 


\dpN/dr 0 ) 


(15b) 


B. Consider Partial Derivatives 

The particular analysis described here includes three 
different consider error sources. These consider parame- 
ters are ( 1 ) the calibration errors in the wet and dry com- 
ponents of the static zenith tropospheric delay; ( 2 ) the cal- 
ibration errors in the daytime and nighttime ionospheric 
delay; and (3) the uncertainties in the station location. All 
delays are expressed in units of length, and fluctuations in 
the delays are ignored here. Although the uncertainties in 
universal time and polar motion (UTPM) are not included 
explicitly, these errors give signatures similar to uncertain- 
ties in station location; for instance, an error in universal 
time is exactly equivalent to an error in station longitude. 
Therefore, the analysis should give a reasonable picture 
of the effects of UTPM errors provided that the assumed 
magnitudes of the station-location errors are large enough 
to encompass the UTPM errors as well. 

The partial derivatives of the Doppler observable with 
respect to calibration errors in the wet and dry components 
of the zenith troposphere have been computed by Russell 2 
and are reproduced here. Define 7 to be the elevation 
angle, and 7 to be its time derivative. Then, for a single 
tracking station, the partial derivative with respect to an 
error Ap Zw in the calibration of the wet zenith delay is 
given by 


df> 

d(Apz w ) 


2/Ma, *2,7,7) 


(16a) 


1 R. K. Russel], “The Development of Parameter Estimation," Engi- 
neering Note 001-100 (internal document), Jet Propulsion Labor- 

atory, Pasadena, California, January 19, 1987. 


2 R. K. Russell, “Computation of Troposphere Partial Derivatives,” 
JPL Engineering Memorandum 391-277 (internal document), Jet 
Propulsion Laboratory, Pasadena, California, February 3, 1972. 
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where 


f(A,B,y,y) = -7 sin 7 + 


tan 7 + B 


x cos 7 — 


(sin 7 -f B cos y) 2 \ 


(16b) 


height of 350 km (z = sin -1 [0.948 cos 7 ]); x is the mini- 
mum of ?r /2 and \2n{t — 0 )/P|, where Pis the period of the 
diurnal variation of the ionosphere (taken to be 32 hours), 
9 is the phase of that variation (taken to be 14 hours of 
time), and t is the local time at the subionospheric point. 
To compute the partials, this time can be written in terms 
of the station location, the universal time, and the direc- 
tion to the spacecraft; the interested reader is referred to 
the work by Jacobson 3 or to [14] for further details. 


The partial with respect to an error A p ZD in the calibration 
of the dry zenith delay is 


dp 

d(Ap ZD ) 




(17) 


Making use of Eq. (5) and the fact that for longitudes 
increasing to the east, d/d\ s = d/d<p (i.e., moving the sta- 
tion to the east is the same as increasing the local sidereal 
time), the partial derivatives with respect to the cylindri- 
cal station coordinates (r,,X s ,z 3 ) easily are shown to be 


where / was defined in Eq. (16b). The constants .4 lt A 2 , 
B 1 , and B 2 were determined by ray-tracing calculations to 
be 


dp 

dr. 


(<p — d) cos 8 sin (y> — a) 
+ 6 sin 6 cos (<p — a) 


Ay = 0.00143' 

A 2 = 0.00035 

> 

By = 0.0445 

B 2 = 0.017 . 


(18) 


dj> 

dX 3 


dp 

dz. 


r, [(<p — d) cos 8 cos ( <p — a) 

— <5 sin 6 sin (ip — a)] 

— 6 cos 8 


> (20) 


The consider partials that relate to errors in calibra- 
tion of the ionospheric propagation have been given by 
Jacobson 3 in terms of the daytime and nighttime “iono- 
spheric coefficients” (V and A f), which essentially repre- 
sent the delay errors at zenith during the day and night. 
The expressions for the partial derivatives are 


In Eq. (20), the partials with respect to r, and z, are 
in units of velocity divided by distance, while the partial 
with respect to A, is in units of velocity per radian. To 
convert the partial with respect to X, to the same units as 
the other partials, the partial derivative given in Eq. (20) 
should be divided by r,. 


= — cos a?— (secz) — x sin x sec z (19a) 

and 

MT = “ 5 (secz) (19b) 

In Eqs. (19a) and (19b), z is the zenith angle of the space- 
craft as seen from the point where the line of sight from 
the tracking station passes through the mean ionospheric 


3 R. A. Jacobson, “Troposphere and Ionosphere Delay Models in 
ATHENA,” JPL Engineering Memorandum 31 J-290 (internal doc- 
ument), Jet Propulsion Laboratory, Pasadena, California, Novem- 
ber 16, 1982. 


V. Weighting Functions 

The consider covariance analysis procedure described 
above carried the weighting factor wy for each Doppler 
measurement; ivy nominally is chosen to be l/Vp 2 in con- 
ventional least-squares estimation. However, Doppler data 
can be affected by large consider errors at low elevations 
because of the effects of the ionosphere and the tropo- 
sphere on signal propagation, as well as the elevation- 
independent effects of uncertain calibration of station lo- 
cations and Earth-orientation parameters. In many cases, 
the inclusion of low-elevation data actually has resulted 
in larger uncertainties in operational navigation than have 
been obtained by discarding those data. Therefore, one 
common approach to this problem has been to discard 
all data below a specified elevation angle. A second ap- 
proach has been to deweight all data uniformly, which 
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has the effect of making the high-elevation data “as bad” 
as the low-elevation data and is not likely to give the 
best trajectory determination. A third possibility, con- 
sidered below, is to deweight the low-elevation data ac- 
cording to some function appropriate to minimize tropo- 
sphere and/or ionosphere effects on the estimated trajec- 
tory. (This deweighting might be accompanied by impo- 
sition of an elevation cutoff.) The single most important 
motivation for constructing the consider analysis method- 
ology described in this article was to make a systematic, 
quantitative investigation of various weighting functions 
that would lead to recommendations regarding appropri- 
ate methods of weighting 8.4-GHz and 32-GHz Doppler 
data. 

The consider analysis described here has been imple- 
mented in a computer program that includes a weighting 
function of the following form: 


Wi 



T > Train 


} ( 21 ) 


= 0 , 


T ^ Tmin 


This weighting function was postulated because it resem- 
bles, in an approximate form, the behavior of the tropo- 
sphere partial derivatives given by Eqs. (16) and (17). As 
before, 7 is the elevation angle of the spacecraft as viewed 
from the tracking station. The elevation-dependent part 
of the weighting function has the form (a t j sin ? 7) 2 , where 
q is a selectable exponent and a e is a selectable weighting 
factor. The assumed Doppler data noise variance is <tx> 2 , 
and all data at elevations below 7 m j n are discarded (this 
elevation cutoff also is user-selectable). Note that with 
this form of weighting, the data analyst has the option 
of adjusting the weighting function in several ways. Data 
can be deweighted uniformly by increasing the value of ao 
and choosing a t = 0. Data can be deweighted as a func- 
tion of elevation angle by appropriate selection of a e and 
q. Finally, 7 min can be varied to change the cutoff angle. 

Figures 4(a) and 4(b) are plots of the partial deriva- 
tives given by Eqs. (16)— (19) as functions of the elevation 
angle 7. These plots are for a segment of the Mars Ob- 
server trajectory, whose state vector is given in Table 1. 
Tracking from the Goldstone Deep Space Network com- 
plex (latitude of about 35.2 deg) is assumed. Note that 
the troposphere partials are dominant at the low elevation 
angles, reaching values about 40 times greater than the 
maximum ionosphere partials in this particular case. For 


observations at 8.4 GHz, the uncertainties in the coeffi- 
cients V and M are comparable to the values of A p ZD and 
Ap Zw , implying that the troposphere calibration error will 
dominate the effects of the ionosphere calibration error in 
the Doppler data. At 2.3 GHz, V and M are larger by a 
factor of ~ 13 and can be much more important. 

Figure 5 shows the Doppler data weights as a function 
of elevation angle for several different values of the variable 
parameters in Eq. (21); of course, the weights at the lower 
elevation points can be set to zero arbitrarily by making a 
choice of 7 m j n . In this figure, an arbitrary value of unity 
is assigned to the Doppler data noise (i.e., <td = 1). Val- 
ues of 1, 2, and 3 are used for q, while cr e is set equal to 
either <td or 0.5 &d ■ For comparison, Fig. 6 is a plot of 
the inverse of the troposphere partial derivatives. If the 
troposphere were the only contributing calibration error 
and dominated the data noise in determining the Doppler 
error, the ideal weighting function would have a shape sim- 
ilar to the curves in Fig. 6. A comparison of Figs. 5 and 6 
implies that the weighting functions with values of q = 2 
or q = 3 may be more appropriate for Doppler data whose 
errors are dominated by errors in the troposphere calibra- 
tion. 

Note that the shapes of the weighting functions in 
Fig. 5 already lead to an interesting conclusion. When the 
weighting function is a strong function of elevation angle 
(i.e., q ss 3), the choice of elevation cutoffs should make 
very little difference in the final results of the analysis. 
This is because, for the larger values of q, the low-elevation 
data already are downweighted so much that they have 
little effect on the orbit determination; they effectively 
have been discarded even before the elevation limit is ap- 
plied. The specification of an elevation cutoff would be 
more important if the ratio of <r e to <tq in Eq. (21) were 
much smaller than the values used in generating the plot 
in Fig. 5. 


VI. Accuracy Over Long Data Arcs 

It was mentioned previously that the trajectory repre- 
sentation developed herein breaks down for long data arcs, 
but the maximum length of time for which the approxima- 
tions are valid was not specified. However, it is possible to 
make a rough estimate of the longest data arcs that might 
be investigated using the formalism described above. This 
estimate has been made by considering the accuracy of the 
linear approximation given by Eqs. (7) and the accuracy of 
the piecewise-linear approximation described by Eqs. (8) 
and (9) and in the surrounding discussion. 
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A more accurate system of equations would include 
terms that are quadratic in time: 


parameters as shown in Table 1. Using the appropriate 
numbers, Eq. (23c) becomes 


r 

6 

a 

r 

6 

a 


— ro + rot + r'ot 2 /2 

— Sq + $ot + Sot 2 j 2 


dr 

dr 0 


1 + 2.6 x 10 _1 V 


(24) 


— oto + dot + dot 2 / 2 
= r 0 + rot + t 2 (d 3 r 0 /dt 3 )/ 2 


} ( 22 ) 


— So 4- Sot 4- t 2 (d 3 S 0 /dt 3 )/2 


= d 0 + d 0 t + t 2 (d 3 a 0 /dt 3 )/2 ) 


where t is given in seconds. Thus, the linear approxima- 
tion used in the standard derivation of the state transition 
matrix A 2 would have a fractional error of 2.6 x 10 -14 < 2 . 
The partial derivative would be accurate at the 1-percent 
level only if t were smaller than 6.2 x 10 5 seconds, or about 
7 days. For 10-percent accuracy in the partial, t could be 
no greater than 23 days. 


As an example, consider the partial derivative dr/dro- In 
the linear theory, this partial derivative is found as follows: 


dr 

dr 0 


1 + tp- = 1 

dr 0 


(23a) 


since ( dr 0 /dr 0 ) = 0. Keeping the quadratic time terms as 
in Eq. (22), letting r g be the line-of-sight acceleration due 
to the Sun’s gravitational field, and using the result from 
[10] that 


ro = r 0 (So + dp cos 2 Sq) 4- r g (23b) 


this partial derivative is 


dr dr 0 t 2 drp 

dr 0 dr 0 2 dr 0 

= 1 + y(^o+«o cos2 ^ 0 + |^) 



} (23c) 


x ( <5o+do cos 2 $0 4 — -~~z [2 — 3 sin 2 ipo 

\ r p,o 

In the above equation, fi is the universal gravitational con- 
stant multiplied by the mass of the Sun, r Pi 0 is the Sun- 
spacecraft distance, and V’o is the Earth-spacecraft-Sun 
angle. The value of dr g /dro has been taken from [10]. 
As an example, a segment of the planned Mars Observer 
trajectory on July 22, 1993 was taken, with epoch state 



The piecewise-linear approximation to the trajectory 
includes the bulk of the effect of the quadratic terms. Con- 
sidering the transition matrix connecting day 2 to day 1, 
Al, the maximum value of t that would be used in Eq. (24) 
would be about a day, or 8.6 x 10 4 seconds. The connec- 
tion from day 3 to day 1 is based on two values of the 
partial derivative, which differ because of the effects of 
the quadratic terms included in a piecewise-linear man- 
ner. Thus, the fractional error for each day’s pass would 
be indicated by Eq. (24), but the total fractional error over 
an n-day data arc is approximately 


/„ « (1 + 2.6 x 10- 14 T 2 )" - 1 

« 2.6 x 10“ 14 nT 2 


(25a) 


instead of 


/„ « [1 + 2.6 x 10 -14 (nT) 2 ] - 1 

= 2.6 x 10~ 14 n 2 r 2 


(25b) 


where T is the length of a day in seconds. These two equa- 
tions differ in only one respect, their dependence on the 
number of days n in the data arc. For n = 7, Eq. (25b) 
predicts a fractional error of 0.95 percent for the linear 
approximation, while Eq. (25a) predicts a fractional er- 
ror of 0.14 percent for the piecewise-linear approximation. 
For a month-long data arc, n = 30, the fractional error for 
the linear approximation would be about 17 percent, while 
that for the piecewise-linear approximation would be only 
0.6 percent! In fact, even this result is an overestimate of 
the inaccuracy of the piecewise-linear approximation used 
in the analysis developed in this article. When the transi- 
tion state matrix between passes j and j — 1 is evaluated 
for the state vector rjj - 1 at the midpoint between the two 
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passes, a more appropriate value of T used in Eq. (25a) 
would be half the length of the day rather than the full 
length of a day. In any case, it is clear that the length of 
the data arc that can be considered is substantially longer 
than that possible for the simple linear approximation. 

Provided that calculations of the neglected terms in 
other partial derivatives give results similar to those de- 
scribed above, for the particular sample trajectory consid- 
ered here, the piecewise-linear approximation represented 
by Eq. (8) should be quite good for data arcs longer than 
a month. The approximation would hold for a longer pe- 
riod cruise trajectory in the outer solar system, say at the 
distance of Jupiter, since the numerical coefficient of T 2 
in Eq. (25a) would be considerably smaller. For a tra- 
jectory near Venus, the gravitational acceleration due to 
the Sun would be larger, reducing the accuracy of the lin- 
ear approximation. In that case, the advantage of the 
piecewise-linear approximation would be even more pro- 
nounced. It is possible to get arbitrarily close to the true 
model of state-vector evolution by breaking the data arc 
into more pieces, but the above estimate shows that the 
simple piecewise-linear approximation that is used should 
be quite good for most circumstances. 

VII. Summary 

This article has described an error covariance analy- 
sis methodology for Doppler measurements that can be 
used to predict errors in the estimated state vector of an 
interplanetary spacecraft in the cruise portion of its tra- 
jectory. The analysis and associated computer program 


improve on previously published analyses in several ways: 
(1) The effects of imperfect calibration of the static tro- 
posphere, the static ionosphere, and the station location 
are all considered simultaneously; (2) the spacecraft state 
transition matrix is modelled more accurately, using a 
piecewise-linear approximation to represent the evolution 
of the spacecraft trajectory; and (3) the Doppler data can 
be weighted according to elevation angle in an attempt to 
minimize the covariance associated with determination of 
the spacecraft trajectory. The first and second “improve- 
ments” or better representations also are available in com- 
prehensive software such as the JPL Orbit Determination 
Program [15], but the limited focus of the analysis method 
described in this article makes it preferable for investiga- 
tion of a large number of possible cases of data weighting 
in a short period of time. 


The analysis method described in this article is useful 
for estimating the errors in trajectory determination for 
data arcs longer than a month in the case of a spacecraft 
more than 1 AU from the Sun and far from the gravita- 
tional influence of any planets. In the case of spacecraft 
in the outer solar system (e.g., Voyager and Galileo), the 
reduced gravitational effect of the Sun and the smaller 
angular velocities of the spacecraft should enable use of 
this formulation for a considerably longer data arc; the 
period of validity needs to be evaluated on a case-by-case 
basis. It is likely that the accumulated effects of nongrav- 
itationa], unmodelled accelerations would lead to greater 
errors than those caused by the piecewise-linear inclusion 
of the quadratic terms due to known gravitational effects. 
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Table 1. Initial state vector for a segment of 
Mars Observer trajectory. 


Parameter 


Value 


Year 

Day of year 
Distance, km 
Right ascension (a), deg 
Declination (5), deg 
Radial velocity, km/sec 
da/dt, deg/sec 
dS/dt, deg/sec 


1993 

203 (22 July) 
3.1573 X 10 8 
169.0252 
5.1308 
11.5770 
6.2454 X 10~ 6 
-2.5866 X 10“ 6 



STATE VECTOR (arbitrary units) 



Fig. 1. Station-spacecraft tracking geometry. 



Fig. 3. Geocentric plane-of-sky Cartesian coordinate system. 



PASS NUMBER 


Fig. 2. The effect of different methods of computing a spacecraft 
state transition matrix. 




I . 1 , Z I . ^ . L_ 

0 20 40 60 80 


ELEVATION ANGLE, deg 

Fig. 4. Partial derivatives of the Doppler error with respect to the 
calibration errors in the troposphere and ionosphere delays plot- 
ted against elevation angle: (a) partials of the Doppler observable 
with respect to wet and dry zenith-troposphere calibration errors, 
and (b) partials of the Doppler observable with respect to daytime 
and nighttime ionosphere coefficients. 



20 40 60 80 

ELEVATION ANGLE, deg 


Fig. 5. The weighting function given by Eq. (21) for a cutoff el- 
evation of 6 deg: Three different values of the exponent q are 
included for each family. 



Fig. 6. Inverse of the partials of the Doppler error with respect 
to errors in the calibration of the zenith troposphere delay. Plots 
for wet and dry components overlap. The plot goes to infinity at 
transit, where -y goes to zero in Eq. (16b). 
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Appendix 

Summary of Partial Derivatives 


This appendix gives (without derivation) the partial 
derivatives in the 6x6 matrices A 2 and A3 described in 
Section III. A2 is the matrix of partials of the state vector 
r at an arbitrary time t with respect to the state vector Fo 
at the epoch time for the same pass. To find the matrix of 
partials relating the state vector at the epoch time on day 
j to that at the epoch time on day j — 1, denoted Aj -1 in 
Section III, F 0 should be replaced by fjj - 1 and t by r, as 
described following Eq. (9). 


4 in 

r Pi 0 = [r 0 2 + r e>0 2 - 2r 0 r ei0 cosxo] ' (A-4) 

The angle <r 0 is given by 

cos (To — — sin 60 cos 6 Si o cos (ao — a,,o) 

+ cos 60 sin 6„ ,o (A-5) 


The matrix A2 is defined as follows: Reference [10] gives the gravitational accelerations of 

the spacecraft in the radial, declination, and right- 
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dr 

dr 

dr 
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ascension directions ( r g ,6g,dg ). If fi is the product of the 

w 
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universal gravitational constant and the mass of the Sun, 

36 
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36 
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those accelerations are 
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3d 
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3d 
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V 3r a 

W 0 

da 0 
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36 0 
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ro 


A number of definitions must be made to write the par- 
tial derivatives in a reasonably compact form; the defini- 
tions made here largely follow those in [10], and all vari- 
ables are assumed to be evaluated at the reference (or 
epoch) time indicated by the subscript “0.” Let the right 
ascension and declination of the Sun be denoted by a,,o 
and 6, 0) respectively. Denote the Earth-Sun distance by 
r C| 0 and the Earth-spacecraft distance by r Pi o- The Sun- 
Earth-spacecraft angle is Xo> where 


x 




(A-6c) 


The nonzero partial derivatives in the matrix A2 are given 
below: 


dr 

dr 0 


= 1 


(A-7a) 


cosxo = cos 60 cos6,,o cos (e*o — «,,o) 

+ sin 60 sin 6, ,o ( A-2) 

The Earth-spacecraft-Sun angle V’o is found from 

sin V'o = *’e,o s i n Xo/ r Pl o (A-3) 


where 



(A-7b) 
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(A-8a) 



(A-8b) 


46 



da 

da 0 

da 
da o 


= 1 


= t 


(A-9a) 
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— = — - — - — 3 — - + (6l + al cos 2 S 0 )t (A-lOa) 
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The matrix A 3 is used to make the transformation be- 
tween a spherical coordinate system and a Cartesian plane- 
of-sky coordinate system. It is defined as follows: 
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The nonzero entries in the matrix A 3 are as follows: 
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